l/l 


AD-A1J1  MC 

UNCLASSIFIED 


EFFICIENT  HIGH  RESOLUTION  CONFUTATIONS  OF  PROJECTILE 
AND  FISSILE  FLON  FIELDS(U)  FEDA  CORF  PALO  ALTO  CA 
C  K  LONAARD  ET  AL.  DEC  *7  AAO-21247.  C-EQ-S 
DAAG29-04-C-MA2  F/O  24/4 


^CCTPTPn 


#■  V/v " 


iC'. 


AD-A191  086  REPORT  DOCUMENTATION  PAGE 

Jib.  RESTRICTIVE  MARKINGS 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION  / DOWNGRADING  SCHEDUL 


a  PERFORMING  ORGANIZATION  REPORT  NUM8ER(S) 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER 

ftfio 


6a.  NAME  OF  PERFORMING  ORGANIZATION 


6b.  OFFICE  SYMBOL  |  7a.  NAME  OF  MONITORING  ORGANIZATION 
(If  applicable)  j 


_ PEDA  Corporation  | 

6c.  ADDRESS  (City,  State,  and  ZIP  Code)  ~ 

4151  Middlefield  Road,  Suite  7 
Palo  Alto,  CA  94303 


U.  S.  Army  Research  Office  ^  )g| 

7b.  ADDRESS  (City.  State,  and  ZIP  Code)  MLtr 

P.  0.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 


8b.  OFFICE  SYMBOL  9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable) 

DAAG29-84-C-0002 


PROJECT 

TASK 

NO 

NO. 

WORK  UNIT 
ACCESSION  NO. 


U.  S.  Army  Research  Office  |  DAAG29-84-L-UUU2 

8c  ADDRESS  (City,  State,  and  ZIP  Code)  10.  SOURCE  OF  FUNDING  NUMBERS 

P,  0.  Box  12211  PROGRAM  PROJECT  j  TASK  IWORK  UNIT 

Research  Triangle  Park,  NC  27709-2211  SimtNT "°  N0  NO-  ACCEWOh  NO 

1 1  TITLE  (include  Security  Classification)  ~~~ 

EFFICIENT,  HIGH  RESOLUTION  COMPUTATIONS  OF  PROJECTILE  AND  MISSILE  FLOW  FIELDS 

12  PERSONAL  AUTHOR(S) 

_ C.K.  Lombard,  Jorge  Bardina,  R.C.-C.  Luh  and  Ethiraj  Venkatapathy 

13a.  TYPE  OF  REPORT  13b.  TIME  COVERED  ^  ,  Jl4.  DATE  OF  REPORT  (Year.  Month,  Day)  IlS.  PAGE  COUNT 

FINAL  from 84  3  15  to  87  9  14  87  12  22  31 

16  SUPPLEMENTARY  NOTATION 

The  view,  opinions  and/or  findings  contained  in  this  report  are  those 
-^l^aA^l-gnnd  Mn?i£g  g^g^Onstl^d  Sfwgg  ^U^h^Ttment  °f  the  Arn,y  P^ition, 

—  7. _ _ COSATI  CODES _ 18.  SUBJECT  TERMS  (Continue  on  reverst  if  necessary  and  identify  by  block  number) 

- field - group - sub-group -  Algebraic  Grid  Generation,  Navier-Stokes ,  Upwind  Methods, 

-  Projectile  Flow,  Missile  Flow,  Base  Flow,  Adaptive  Grid, 

- - - - Composite  Grid.  Overset  Grid _ 

'9  ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

J  A  new  family  of  computational  tools  is  explored  in  the  context  of  generic  problems  of  interest  to  Army 
Research  and  Development  missions.  The  problems  are  the  SOCBT  projectile  flow  and  the  MICOM  model 
tactical  base  flow.  Both  problems,  for  which  experimental  data  is  available  for  comparison,  have  also  been 
extensively  computed  by  other  groups  with  other  methods,  principally  the  central  difference  method  of  Pp«m 
and  Warming. 

Here  the  techniques  we  employ  and  to  some  extent  refine  are  composite  and  overset  patched  mesh 
systems  generated  with  a  new  multistep,  graphics  interactive  algebraic  procedure  FASTWO  and  the  CSCM- 
S  symmetric  Gauss-Seidel  relaxation  procedure.  The  grid  generation  technique  permits  comparatively  simple 
direct  user  control  of  mesh  spacing  and  quality. 

The  CSCM-S  algorithm  is  an  alternating  direction  space  marching  method  of  lines  in  2-D  and  planes 

IQ  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

□  unclassified/unlimited  □  SAME  AS  RPT,  □  otic  USERS _ Unclassified _ : 

22a  NAME  OF  RESPONSIBLE  INDIVIDUAL  22b  TELEPHON’  (Include  Area  Code)  |  22c  OFFICE  SYMBOL 


OO  FORM  1473, 84  mar 


83  APR  edition  may  be  used  until  exhausted 
All  other  editions  are  obsolete 


SECURITY  CLASSIFICATION  OF  this  page 
UNCLASSIFIED 


"  'TUHl^L*  ILFltFlLF  V*  WK.^  X.”  V  7~.V  kwjva'v.  irwi  ■V-  7vwwi«'..-w  -. 


vSSIFIED 


;  .  CLASSIFICATION  OF  THIS  FAOt 


V 


19.  ABSTRACT  continued 


in  3-D.  The  method  is  economical  both  in  storage  and  data  handling  and,  also,  is  rapidly  convergent  as  the 
result  of  improved  propagation  of  nonlinear  advection,  boundary  to  boundary.  In  the  SOCBT  problem  the 
method  exhibited  convergence  in  order  100  global  iterations  and  showed  improved  accuracy  vis  a  vis  central 
difference  methods  in  capturing  features  of  the  flow  such  as  shocks  and  expansions.  The  method  is  stable  in 
both  elliptic  regions  and  zones  of  streamwise  separated  flow  including  base  flow. 
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Figure  Titles 


Figure  1  SOCBT  Grid  Generation 

a.  Patch  Boundaries 

b.  Complete  Grid 

c.  Near  Field 

d.  Base  Region 

Figure  2  Mach  Contour  Plots  for  SOCBT  at  Mach  1.1  Freestream 

a.  Boattail  Continued  Afterbody 

b.  Sting  Mounted  Base 

Figure  3  Cp  Plots  for  SOCBT  at  Mach  1.1 

a.  Boattail  Continued  Afterbody 

b.  Sting  Mounted  Base 

Figure  4  Plots  for  SOCBT  at  Mach  0.98 

a.  Mach  Contours  Free  Flight  Base 

b.  Cp,  Boattail  Continued  Afterbody 

c.  Cp,  Sting  Mounted  Base 

d.  Cp,  Free  Flight  Base 

Figure  5  Perspective  View  of  the  SOCBT  Surface  Grid  for  the  Symmetry  Half  Body 

Figure  6  Cp  Plots  for  SOCBT  at  Mach  0.96,  Angle  of  Attack  4° 

a.  CSCM-S 

b.  Sahu51 

Figure  7  Mach  Contour  Plot  SOCBT  at  Mach  0.96,  Angle  of  Attack  4° 

Figure  8  Cp  Plots  for  SOCBT  at  Mach  0.95,  Angle  of  Attack  2° 

a.  CSCM-S 

b.  Nietubicz,  et  al52 

Figure  9  Grids  for  MICOM  Base  Flow  Study 

a.  Composite,  Double  Wraparound 

b.  Composite,  Cartesian  Step 

c.  Overset,  Multiple  Patch 

d.  Solution  Adaptive  from  (c) 

Figure  10  Mach  Contour  Plots  for  MICOM  Base  Flow  Computed  with  CSCM-S 

a.  Composite,  Double  Wraparound 

b.  Composite,  Cartesian  Step 

c.  Overset,  Multiple  Patch 

d.  Solution  Adaptive  from  (c) 

Figure  11  Base  Pressure  Plots  for  MICOM  Base  Flow  Computed  with  CSCM-S 

a.  Composite  Grid  Results 

b.  Overset  Grid  Results 

Figure  12  Base  Pressure  Plots  for  MICOM  Base  Flow  from  Earlier  Studies15  with  Other 
Methods 
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1.  Introduction 


Based  on  the  Beam- Warming  central  difference  algorithm1,2  of  the  late  sev¬ 
enties  and  with  the  increasing  power  of  scientific  computers,  the  early  eighties 
saw  a  substantial  computational  effort3- 10  at  analyzing  projectile  and  missile  flow- 
fields.  The  early  work3-9  focused  particular  attention  on  the  important  Secant- 
Ogive-Cylinder-Boattail  (SOCBT)  projectile  configuration  for  which  a  substan¬ 
tial  body5'6,7,17  of  experimental  results  had  been  obtained.  The  results  of  an 
assessment18  made  shortly  before  the  present  work  was  initiated  may  be  summa¬ 
rized  as  follows.  At  zero  angle  of  attack  or  away  from  the  immediate  transonic  zone 
.9  <  M  <  1.1,  the  subsonic  and  supersonic  techniques  using  the  thin  layer  Navier- 
Stokes  approximation19  in  respective  time1  and  space  marching2  implementations 
of  Beam  and  Warming  factored  implicit  numerical  methods  had  been  found  effective 
in  flow  field  predictions  and  the  derived  aerodynamic  coefficients.  In  the  immediate 
transonic  regime  where  aerodynamic  coefficients  are  a  highly  nonlinear  and  rapidly 
varying  function  of  Mach  number,  performance  had  not  been  as  satisfactory.  In 
the  later  regime,  the  variation  of  coefficients  is  largely  governed5,7  by  the  sensitive 
movements  of  shocks,  aftward  with  increasing  Mach  number;  and  computational 
accuracy  and  numerical  stability  become  of  major  importance.  Intimately  related 
to  the  shocks  are  precursor  expansions  which  are  equally  difficult  and  important  to 
treat  accurately.  A  final  computational  detail  needing  improvement  was  the  ability 
to  accurately  compute  flow  in  the  vicinity  of  contact  discontinuities  which  appear 
in  streamwise  and  cross  flow  separations,  for  example. 

Another  emerging  computation  area  was  base  flow  analysis  which  became  an 
important  research  topic  for  both18  the  impacts  on  drag  and  pitching  moment.  Base 
flow  computations  were  performed  for  projectiles7  and,  with  more  emphasis10-16, 
on  a  wind  tunnel  model  tactical  missile  with  centered  propulsive  jet.  The  latter 
computations  provided  universally  unsatisfactory  predictions  of  base  pressure,  a 
deficiency  largely  ascribed  at  the  time  to  inadequate  turbulence  modeling. 

In  the  same  time  frame  of  these  early  computations,  PEDA  Corporation  was 
carrying  out  research  and  development  to  provide  advanced  computational  tools 
for  viscous  flow  analyses  in  complicated  geometries.  The  theme  of  the  research 
was  accurate,  efficient  flux  difference  split  implicit  upwind  schemes20-25  that  did 
not  require  ad  hoc  smoothing  operators  and  efficient  high  quality  algebraic  grid 
generation26  on  domain  boundary  fitting  patch  mesh  systems. 

The  flux  difference  split  upwind  implicit  schemes  with  coupled  well  posed  char¬ 
acteristic  boundary  point  approximations  were  known  to  have  superior  natural  sta¬ 
bility  properties  while  providing  sharper  capture  of  flow  structures  such  as  shocks, 
expansions  and  contacts.  A  newly  discovered  symmetric  Gauss-Seidel  relaxation 
scheme  CSCM-S24,25  for  the  CSCM  upwind  method20-93  had  been  found  to  be 
both  more  efficient  in  data  storage  and  handling  but,  also,  to  be  substantially  more 
rapidly  convergent  to  steady  state. 
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The  complementary  development,  algebraic  grid  generation  in  patched  mesh 
systems,  had  the  fundamental  capability  to  provide  direct  control26,27  of  mesh  topol¬ 
ogy,  spacing  and  quality  without  the  necessity  and  numerical  complexity  of  solving 
sets  of  partial  differential  equations,  as  was  being  practiced28  in  some  prominent 
lines  of  development. 


The  intent  of  the  presently  reported  effort  was  to  implement,  test  and  demon¬ 
strate  the  applicability  and  effectiveness  of  the  evolving  new  computational  method¬ 
ology  in  the  analysis  of  problems  of  interest  to  US  Army  research  and  development. 
Based  on  both  theory  and  numerical  experience,  the  belief  was  that  the  new  methods 
had  the  potential  to  improve  accuracy,  efficiency  and  productivity  in  aerodynamic 
design  analysis. 

In  the  next  section  we  will  summarize  the  substantial  effort  that  has  been  made 
in  adapting  tools  and  exploring  the  numerical  solution  of  two  generic  transonic 
aerodynamics  problems  -  SOCBT  and  MICOM  missile  base  flow  -  of  direct  Army 
interest  and  for  which  good  experimental  and  computational  background  exist. 

2.  Research  Accomplished 

Over  the  nominally  three  year  program  we  have  explored  the  2-D  axisymmetric 
problems  of  the  SOCBT  forebody  and  base  flow  and  MICOM  missile  base  flow. 
That  work  has  resulted  in  sue  meeting  papers29"34  reporting  different  advances  and 
understandings  obtained  in  the  course  of  attacking  these  problems.  We  have  also 
studied  3-D  flows  about  the  SOCBT  at  angle  of  attack.  In  the  latter  part  of  the 
program  we  requested  and  obtained  a  six  month  no  cost  extension  permitting  us 
under  other  support  to  further  advance  our  3-D  computational  tools  -  algebraic 
grid  generation  and  the  CSCM-S  algorithm.  The  latter  improved  program35  has 
been  employed  in  the  3-D  angle  of  attack  computations  for  the  SOCBT  that  we 
report  on  here  and  elsewhere36. 

Algebraic  Grid  Generation 

The  concept  of  patched  meshing  in  which  complex  domains  are  broken  up 
into  many  geometrically  regular  and  topologically  rectangular  subdomains  leads 
naturally  to  the  use  of  efficient  algebraic  techniques  for  the  construction  of  the  in¬ 
dividual  mesh  patches.  To  obtain  the  desired  smoothness  properties  over  the  global 
mesh  in  the  vicinity  of  patch  boundaries,  comparatively  early  in  our  methods  re¬ 
search,  a  technique  that  permits  specification  of  point  distribution  and  gradient  on 
all  boundaries  had  been  devised.  The  technique  -  termed  generalized  transfinite 
interpolation26  -  can  be  seen  to  embed  use  of  a  parameterized  general  cubic  poly¬ 
nomial  for  the  coordinate  curves.  Regularity  of  the  mesh  is  obtained  by  employing 
continuous  distributions  of  the  parameters  of  the  curves  within  judiciously  chosen 
bounds  based  on  analysis.  Stretching  functions  such  as  that  of  Vinokur37  are  used 
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to  distribute  points  and  blending  funct  ms  axe  used  to  distribute  parameters  of  the 
curves  between  lateral  boundaries. 

A  novel  feature  of  the  technique  was  the  introduction  of  the  corner  singularity 
from  analysis  to  govern  distribution  of  points  and  parameters  in  the  vicinity  of 
boundary  slope  singularities.  At  such  points,  the  method  thus  obtains  the  desired 
properties  of  mesh  smoothness  to  the  interior. 

However,  in  attempting  to  apply  the  generalized  transfinite  interpolation  tech¬ 
nique  in  a  variety  of  2-D  problems,  it  became  evident  the  method  was  too  sensitive 
to  parameter  selection  among  too  many  options,  was  confusing  and  ultimately  re¬ 
quired  too  much  artistry  to  meet  the  objectives  of  simplicity  and  user  friendliness 
set  for  the  products  of  the  research.  Further,  the  lack  of  an  analytical  solution  to 
corner  problems  blocked  the  straightforward  extention  of  the  technique  to  3-D. 

With  some  reflection  it  became  evident  to  us  early  in  the  present  program 
that  the  difficulty  lay  in  trying  to  accomplish  too  much  in  a  single  step  process. 
Rather,  borrowing  the  tools  of  the  algebraic  technique  and  redefining  the  process 
in  a  sequence  of  simple  steps  with  interactive  computer  graphics,  we  could  define  a 
straightforward  procedure  to  meet  the  desired  ends. 

The  fundamental  approach  that  we  employ  and  which  was  initially  implemented 
under  other  support  is  in  the  realm  of  two  boundary  methods,  in  that  one  pair  of 
opposite  sides  of  a  patch  is  regarded  as  prescribed  and  often  includes  a  portion  of 
a  physical  boundary.  The  other  pair  of  sides  is  formed  of  the  left  and  right  limiting 
members  of  the  family  of  generalized  cubic  coordinate  curves  joining  the  initially 
given  two  boundaries.  In  either  2-D  or  3-D  the  general  cubic  coordinate  curve  has 
the  simple  form 

L  =  ri  +  (r2-  ri)/(u)  +  0\g{u)  +  a2.h(u)  (l) 

where 

/(u)  =  u2(3  -  2u) 
ff(u)  =  u(l  -u)2 
h(u)  =  u2(u  -  1) 

Equation  (l)  is  a  hermite  interpolation  of  value  (r)  and  gradient  (a)  on  the  two 
boundaries  and  is  parameterized  in  terms  of  u  which  varies  from  zero  to  unity.  The 
scalings  of  and  £2  influence  the  shape  (curvature)  of  the  curve  between  any  pair 
of  end  points.  The  specification  of  a  discrete  set  of  u  values  using  a  generalized 
distribution  function  such  as  that  of  reference  37  defines  the  nodal  intersections 
with  the  other  family  of  coordinate  curves. 

In  our  implementation,  the  left  and  right  limiting  (lateral  bounding)  coordinate 
curves  are  developed  interactively  on  a  graphics  terminal  or  workstation  to  have 
the  desired  configurations.  The  parameters  of  these  lateral  bounding  curves  are 
then  blended  with  polynomial  weighting  functions  to  describe  the  general  cubic 
coordinate  curve  over  the  intervening,  also  discretized,  interval. 
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The  lateral  patch  boundaries  are  essentially  control  devices  that  specify  shape 
and  distribution  to  surrounding  regions.  As  such  they  are  placed  where  needed  - 
at  breaks  in  body  surface  geometry  and  as  terminators  or  transition  guides  from 
regions  of  strong  shape  variation  to  regions  of  very  regular  mesh.  In  fact  the  mesh 
generation  problem,  particularly  for  geometries  with  any  substantial  complexity,  is 
a  problem  of  multiple  length  scales.  The  purpose  of  multiple  patching  is  to  isolate 
regions  of  comparable  scales  and  on  which  subdomains  the  solution  is  comparatively 
regular  and  can  be  conveniently  fit  by  simple  functions. 

Once  a  primary  grid  is  generated  by  the  technique  described  above  it  can  be 
interactively  improved  by  modifying  parameter  blending  and  point  distribution  in¬ 
cluding  point  redistribution  along  the  alternate  family  of  coordinate  lines  implicitly 
defined  by  the  nodes  on  the  cubic  coordinate  curves.  The  latter  operation  is  in 
the  spirit,  if  not  the  detailed  implementation,  of  a  two  step  generalized  transfinite 
interpolation. 

Another  secondary  operation  that  we  employ  is  the  modification  of  coordinate 
lines  in  the  vicinity  of  a  boundary  to  smoothly  enforce  local  normality.  The  latter 
operation  like  all  the  procedures  has  been  programmed  as  a  convenient  tool  requiring 
minimal  input  to  apply  at  a  boundary.  Finally,  a  parameterized  tension  spline27 
that  provides  an  analytical  description  of  a  curve  amongst  discrete  data  is  a  tool  that 
has  proved  useful  in  the  latter  operation,  in  effecting  redistribution  of  points  along 
coordinate  curves  of  either  family  and  for  fitting  numerically  specified  boundary 
data. 

Under  ARO  support  we  have  explored  the  application  of  the  simplified  mul¬ 
tistep  technique  to  the  projectile  and  missile  problems  we  have  attacked  computa¬ 
tionally.  In  the  course  of  the  work  and  particularly  over  the  first  two  years,  we  have 
progressively  identified  ways  of  improving  the  organization  of  the  tools  to  provide 
better  user  control  of  mesh  quality  and  enhanced  productivity.  The  tools  and  orga¬ 
nization  of  patched  grid  generation  in  five  algebraic  steps  are  described  in  references 
33  and  34.  The  five  step  procedure  (FASTWO  in  2-D)  that  has  evolved  is 
Step  1.  Boundary  definition  for  the  global  computational  domain 
Step  2.  Wireframe  block  (patch)  decomposition  of  the  computational  domain 
Step  3.  Trial  grid  generation 
Step  4.  Grid  distribution  adjustments 
Step  5.  Grid  normalization  at  boundaries. 

Over  the  past  three  years,  the  work  has  been  very  effectively  employed  in  other 
Army  R  &  D  programs  -  notably  grids  for  coupled  external/cavity  flow  simula¬ 
tions  in  the  transonic  (AOA  aircraft)39  and  hypersonic  (HEDI  missile)40,41  speed 
regimes  for  USA-SDC.  The  work  is  presently  finding  application  with  USA-BRL  in 
an  analysis  for  sabot  separation. 

Many  examples  of  grids  generated  with  the  technique  are  given  in  the  references 
to  our  work.  Here,  for  illustration,  we  show  in  Figures  la,  lb,  lc  and  Id  the 
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wire  frame  domain  decomposition  for  the  SOCBT  problem  and  three  views  of  the 
resulting  grid. 

SOCBT  Flows  with  the  CSCM-S  Algorithm 

The  technique  that  we  have  applied  under  the  contract  is  a  new  single  level 
operationally  explicit  but  effectively  implicit  algorithm24  for  gasdynamics.  The  al¬ 
gorithm  is  particularly  appropriate  for  multiple  patch  mesh  systems  because  through 
the  interior  boundary  treatment42  each  solution  sweep  operation  on  any  patch  can 
be  decoupled  from  any  other.  Thus  the  method  is  not  only  very  storage  efficient 
and  simple  to  program  including  the  coupling  at  patch  boundaries;  but,  also,  can 
make  excellent  use  of  parallel  computing  in  several  straightforward  ways. 

Previously  the  Beam- Warming  factored  implicit  algorithm1  with  the  Baldwin- 
Lomax  thin  layer  viscous  approximation19  has  provided  the  basis  for  two  similar 
space  marching  (PNS)  procedures2,43  for  the  compressible  Navier-Stokes  equations. 
These  PNS  methods  which  are  highly  efficient  -  requiring  half  the  data  storage  and 
a  small  fraction  of  the  computer  time  of  two  level  time  dependent  methods  -  have 
proven  effective  for  flowslu  with  favorable  streamwise  pressure  gradient  or  with  rela¬ 
tively  small  adverse  pressure  gradients.  However,  in  the  presence  of  a  strong  adverse 
pressure  gradient  such  as  occurs  in  a  wing  or  fin  root  regions  the  contemporary  PNS 
methods  suffer  numerical  stability  problems  and  may  infer  streamwise  separation 
even  where  separation  doesn’t  occur44.  In  such  unseparated  (perhaps  weakly  sep¬ 
arated)  regions,  numerical  stability  has  been  maintained  at  the  price  of  employing 
large  amounts  of  artificial  viscosity  with  a  resulting  loss  in  predictive  accuracy  and 
knowledge  of  the  actual  state  of  the  flow.  Where  strong  streamwise  separation  oc¬ 
curs  the  methods  are  unstable  and  cannot  proceed.  Streamwise  separation  becomes 
a  likely  occurrence  in  fin  roots  and  on  deflected  control  surfaces  and  base  regions. 
Thus  a  more  general  technique  is  needed  that  is  inherently  stable  for  all  types  of 
upstream  influence.  At  a  minimum  the  mixed  elliptic-hyperbolic  problem  requires 
global  iteration,  preferably  with  type  dependent  differencing.  More  background  on 
this  problem  area  is  given  in  reference  24. 

Universal  Single  Level  Scheme  CSCM-S 

The  CSCM  flux  difference  eigenvector  split  upwind  implicit  method20-23  for 
the  inviscid  terms  of  the  compressible  Navier-Stokes  equations  provides  the  natural 
basis  for  an  unconditionally  stable  space  marching  technique  through  regions  of 
subsonic  and  streamwise  separated  flow.  In  such  regions  the  split  method  can  be 
likened  to  stable  marching  of  each  scalar  characteristic  wave  system  in  the  direction 
of  its  associated  eigenvalue  (simple  wave  velocity).  In  supersonic  flow,  where  all 
eigenvalues  have  the  same  sign,  the  method  automatically  becomes  equivalent  to  the 
referenced  PNS  techniques  based  on  the  Beam- Warming  factored  implicit  method 
with  the  Baldwin-Lomax  thin  layer  viscous  approximation. 

Compared  to  contemporary  central  difference  methods,  the  CSCM  character¬ 
istics  based  upwind  difference  approximation  with  its  inherent  numerical  stability 


leads  to  greatly  reduced  oscillation  and  greater  accuracy  in  the  presence  of  captured 
discontinuities  such  as  shocks,  contacts  and  physical  or  computational  boundaries. 
The  method  has  well  validated  heat  transfer  capability45"47.  The  correct  math¬ 
ematical  domains  of  dependence  that  correspond  with  physical  directions  of  wave 
propagation  are  coupled  with  well  posed  characteristic  boundary  approximations23 
naturally  consistent  with  the  interior  point  scheme.  The  result  is  a  faster  sorting  out 
of  transient  disturbances  and  substantially  more  rapid  convergence  to  the  steady 
state.  The  splitting  and  the  associated  time  dependent  implicit  method  have  been 
described  in  detail  in  references  (20)  and  (23)  respectively  for  quasi  1-D  and  2-D 
planar  or  axisymmetric  flow. 

In  the  following,  we  will  sketch  the  differences  between  the  time  dependent 
method  and  the  space  marching  technique  which  we  have  designated  CSCM-S.  The 
discussion  will  begin  with  the  quasi  1-D  inviscid  formulation,  then  give  additional 
details  entering  into  multidimensional  inviscid  and  thin  layer  viscous  procedures 
and  present  some  2-D  axisymmetric  solutions  for  the  SOCBT.  Lastly  we  sketch  a 
3-D  implicit  method  of  planes  algorithm  and  give  some  results  for  the  SOCBT  at 
angle  of  attack. 

Quasi  1-D  Formulation 

The  general  jth  interior  point  difference  equations  for  the  time  dependent 
CSCM  upwind  implicit  method  for  the  inviscid  advection  terms  is 

(7  +  I+V  +  I“A)%  =  ~A+Aq)  n  - A~Aq)n  (2) 

3  ~  1  3 

where  V  and  A  are  backward  and  forward  spatial  difference  operators.  Here  q  is 
the  conservative  dependent  variable  vector  and  F  is  the  associated  flux  vector.  In 
the  notation  the  interval  averaged  matrices  between  node  points  j  and  j  +  1  are 
indexed  j.  For  simplicity,  the  right  hand  side  of  equation  (2)  is  written  for  the  first 
order  method.  Higher  order  methods  in  space  are  given  with  results  in  references 
20,  23  and  38.  In  equation  (2)  the  CSCM  flux  difference  splitting  is 

(A+  +  A~)Aq  =  AF+  +  AF-  =  AF  (3a) 

with 

A±  =  (mTI±T~1'M~1)A  =  A±A  (3  b) 

and 

/±  =  ^(7  ±  sgn(A))  (3c) 

exhibiting  the  similarity  transformation  that  diagonalizes  the  constructed20,23  flux 
difference  Jacobian  A.  Here  A  is  a  diagonal  matrix  of  the  interval  averaged  eigen¬ 
values  that  through  the  truth  function  diagonal  matrices  I ±  make  the  decisions 
about  directions  of  characteristic  wave  propagation  and  whether  or  not  to  send 
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signal  to  the  solution  point.  Thus  in  equation  (2)  the  piece  of  the  flux  difference 
splitting  A+Aq)j-i  represents  the  convection  of  characteristic  wave  contributions 
in  the  positive  coordinate  direction  from  grid  point  j  —  1  to  solution  point  j  and  A~ , 
in  the  negative  direction  from  j  +  1  to  j.  As  the  result  of  incorporating  multiplica- 
tively  the  (local)  time  step  (for  pseudo  time  relaxation)  and  the  spatial  (divided) 
differences  in  the  matrices,  the  numerical  eigenvalues  are  Courant  numbers  for  the 
characteristic  waves  whose  speeds  are  u,  u  +  c  and  u  -  c,  with  c  the  sound  speed. 

Central  to  its  accurate  shock  capturing  capability,  the  CSCM  conservative  flux 
difference  splitting  has  the  Roe48  “property  U”  embodied  in  equation  3a. 

With  6q  =  qn+1  —  qn,  equation  (2)  defines  a  two  level  linearized  coupled  block 
matrix  implicit  scheme  that  can  be  solved  by  a  block  tridiagonal  procedure.  In 
reference  (23)  a  new  (DDADI)  approximately  factored  alternating  sweep  bidiagonal 
solution  procedure  for  equation  (2)  is  presented  that  is  shown  to  be  very  robust 
and  is  operationally  explicit,  i.e.  requires  only  a  decoupled  sequence  of  local  block 
matrix  inversions  rather  than  the  solution  of  the  coupled  set.  For  the  forward  sweep 
the  bidiagonal  solution  procedure  can  be  written 

{I  +  A+  -A~)6q*}.  =  RHS  +  A+Vy-i  (4) 

For  the  linear  problem,  i.e.  constant  coefficient  case  of  stability  analysis,  equation 
(4)  is  equivalent  to  the  single  level  space  marching  procedure 

(/+  A*  -  A~)6q'i  =  A V,-i  -  A*g"  -  A~ A,)"  (5) 

Nonlinearity  enters  in  the  single  level  space  marching  form  (5)  in  that  at  each  step  of 
the  forward  sweep  the  matrices  A+  are  averaged  between  q* y_j  and  q™.  rather  than 
homogeneously  at  the  old  iteration  level  n.  Similarly,  a  companion  backward  space 
marching  sweep  that  is  symmetric  to  equation  (5)  and  that  is  intimately  related  to 
the  backward  sweep  of  the  alternating  bidiagonal  algorithm  of  reference  (23)  is 

(I  +  A+-  1-)%  =  -A+ Aq*)j-i  +  A~q*  -  A~q"  +  J  (6) 

The  method  given  by  equations  (5)  and  (6)  is  von  Neumann  unconditionally  sta¬ 
ble  for  the  scalar  wave  equation.  The  analysis  shows  the  significance  of  DDADI 
approximate  factorization  in  rendering  both  the  forward  and  backward  sweeps  sep¬ 
arately  stable  regardless  of  eigenvalue  sign.  Consequently  as  the  local  Courant 
number  becomes  very  large,  the  robust  method  becomes  a  very  effective  (symmetric 
Gauss-Seidel)  relaxation  scheme  for  the  steady  equations,  a  fact  which  substantially 
contributes  to  the  very  fast  performance  that  will  be  demonstrated. 

At  a  right  computational  boundary  on  the  forward  sweep  we  solve  the  charac¬ 
teristic  boundary  point  approximation21,23 

(i++i+)«,*  (7) 


(8) 


=  q  and  at  a  left,  on  the  backward  sweep 

~  x  „  ~  n  n  +  1 

( A  -A  )6q1  =  A  q  x  —  A  q  ^ 

Following  the  solution  of  equations  (7)  and  (8)  the  conservative  state  vector  is  iter¬ 
atively  corrected21,23  to  maintain  the  accuracy  of  prescribed  boundary  conditions 
while  not  disrupting  the  representation  of  the  computed  characteristic  variables  run¬ 
ning  to  the  boundary  from  the  interior.  Analysis  of  a  model  system  with  upwind 
differenced  scalar  equations  and  coupled  boundary  conditions  was  related  to  the 
linearized  bidiagonal  scheme23  by  Oliger  and  Lombard49;  the  analysis  also  strongly 
supports  the  numerically  confirmed  robust  stability  of  the  present  nonlinear  method 
for  gasdynamics.  A  useful  result  of  reference  49  that  simplifies  the  procedure  of  ref¬ 
erence  23  is  that  on  the  forward  sweep  there  is  no  need  for  a  predictor  step  at  the 
left  boundary  J  =  1;  thus,  the  solution  sweep  begins  at  J  —  2.  Similarly,  the 
backward  sweep  begins  at  J  =  N  —  1. 

With  the  updating  at  each  step,  where  in  equation  (6)  Sqj  =  g”+1  —  it  is 
clear  that  the  symmetric  pair  of  equations  (5)  and  (6)  serve  to  advance  the  solution 
two  pseudo  time  (iteration)  levels;  whereas,  the  linear  alternating  bidiagonal  sweep 
algorithm  of  reference  (23)  advances  the  solution  only  one  level.  To  maintain  con¬ 
servation  to  a  very  high  degree,  in  single  sweep  marching  in  supersonic  zones  we 
iterate  (at  least)  once  locally  at  each  space  marching  step.  The  local  iteration  serves 
to  make  the  eigenvectors  in  the  coefficient  matrices  consistent  with  the  advanced 
state  and  thus  provides  improved  accuracy  for  the  nonlinear  system.  It  appears 
effective  to  do  this  inner  iteration  everywhere,  i.e.  in  both  subsonic  and  supersonic 
regions,  as  the  number  of  global  iteration  steps  to  convergence  with  two  inner  iter¬ 
ations  has  been  found  reduced  by  a  factor  of  three  to  four.  Since  the  computational 
work  per  two  steps  is  about  the  same  for  the  single  level  and  two  level  schemes,  and 
beyond  the  fact  that  one  saves  a  level  of  storage  in  the  space  marching  algorithm, 
the  question  arises:  Can  one  get  solutions  in  less  computational  work  through  faster 
convergence  with  the  nonlinear  space  marching  algorithm? 

With  the  globally  iterative  nonlinear  space  marching  formulation,  early 
experience24  in  two  quasi  1-D  nozzle  problems  with  mixed  supersonic-subsonic  zones 
is  that  solutions  are  obtained  in  roughly  an  order  of  magnitude  fewer  iteration  steps 
than  had  been  required  with  the  previously  fast  pseudo  time  dependent  technique 
and  block  tridiagonal  solving. 

Two  Dimensional  Formulation 

For  two  dimensional  flow,  assuming  a  marching  coordinate  £,  inviscid  terms 

B+Vn  +  B~  A„  (9a) 

and 

-B  +  A^)*_!  -B~Anq)k  (96) 
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are  added  to  the  left  and  right  hand  sides  respectively  of  both  the  forward  and  back-  J 

ward  sweep  equations  (5)  and  (6).  For  viscous  flow,  second  centrally  differenced,  j 

thin  layer  viscous  terms  are  also  added  in  the  r?  direction  as  is  conventionally  prac¬ 
ticed,  e.g.  Steger50.  With  the  terms  for  the  rj  cross  marching  coordinate  direction, 
the  technique  now  becomes  an  implicit  method  of  lines.  Along  each  r]  coordinate 
line,  one  can  solve  the  equations  coupled  with  a  block  tridiagonal  procedure.  Alter¬ 
natively,  a  further  DDADI  bidiagonal  approximate  factorization  can  be  employed  in  j 

the  r)  direction  and  solved  either  linearly  as  in  reference  (23)  or  nonlinearly  as  here 
in  the  £  direction.  As  shown  in  the  quasi  2-D  numerical  experiments  of  reference  I 

(23),  DDADI  bidiagonal  approximate  factorization  is  stable  for  viscous  as  well  as 
inviscid  terms.  Finally  in  reference  (23)  there  is  a  relevant  discussion  of  the  reduced 
approximate  factorization  error  that  attends  using  DDADI  in  one  or  more  space 
directions. 

SOCBT  Axisymmetric  Flow 

In  the  course  cf  the  computational  analyses29  we  simulated  viscous  transonic 
flows  at  Mach  numbers  .9,  .94,  .96,  .98  and  1.1.  Various  approximations29  of 
base/wake  treatment  with  respect  to  the  free  flying  projectile  and  wind  tunnel 
sting  mount  were  undertaken.  These  can  be  seen  in  the  Mach  contour  plots  of  Fig¬ 
ures  2a  and  2b.  In  supersonic  flow,  for  the  case  of  the  continued  boattail  sans  base 
flow/wake  region,  Figure  2a,  one  loses  the  detail  of  a  corner  expansion  followed  by 
a  lip  shock  that  occurs  with  a  backward  step,  Figure  2b.  For  comparison,  Cp  plots 
for  the  two  cases  are  shown  in  Figures  3a  and  3b.  The  low  pressure  of  the  base 
expansion  is  lost  with  the  extended  sting  approximation. 

For  the  high  subsonic  regime,  Figures  4a,  4b,  4c  and  4d,  the  base  pressures 
can  be  seen  to  be  slightly  underpredicted  for  even  the  sting  mounted  case,  Figure 
4c,  vis  a  vis  the  free  flight  geometry,  Figure  4d.  In  both  supersonic  and  subsonic 
flows,  Figures  3a  and  4b,  the  upwind  differencing  scheme  produced  qualitatively 
sharper  responses  to  flow  structures  and  in  better  agreement  with  experiment  than 
the  earlier  computations  with  central  differencing,  with  the  attendant  smoothing. 

More  details  concerning  the  experiment  and  computational  results  can  be  found 
in  references  17,  5  and  29. 

Three  Dimensional  Method  of  Planes  Algorithm 

In  reference  24  we  presented  a  symbolic  algebra  for  DDADI  approximate  fac¬ 
torization  and  derived  single  level  relaxation  schemes.  The  algebra  is  based  on  the 
implicit  difference  stencil  of  the  implicit  method.  Here  we  use  the  approach  to 
derive  the  symmetric  Gauss-Seidel  implicit  method  of  planes  relaxation  algorithm 
employed  in  the  present  work. 

The  unfactored  three  dimensional  linearized  implicit  method  can  be  represented 
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by  the  symbolic  matrix  expression 


-C+ 


B~ 

C~ 

D  A' 

-B+ 


On  the  block  diagonal  the  matrix  D  =  I  +  A+  -  A~  +  B+  —  B~  +  C+  —  C~ .  A 
once  DDADI  approximate  factorization  in  the  r\  coordinate  direction  leads  to  the 
expression 


•  D- 


-C+ 


—B+ 


-C+ 


—B+ 


By  analogy  with  the  derivation  of  the  single  level  scheme  for  quasi  1-D  flow  from 
DDADI  bidiagonal  approximate  factorization,  we  identify  the  above  expression  with 
the  alternate  space  marching  implicit  method  of  planes  algorithm 

Forward  Sweep 


-c+  D  C-  ] 

-B*  J 

Backward  Sweep 


-c+  D  c-  Sg]*1  =  ,;+1 , ] 

-B+ 


In  the  planes  the  coupled  block  matrix  problem  can  be  further  simplified  by  the 
approximate  factorization 


[  -C+  ,  D  ,  C~  )  D~l  [  -B+  D  B~  ]6q  =  RHS 
which  leads  directly  to  the  block  tridiagonal  solution  sequence 

[  -C+  D  C~  }Sq *  =  RHS 
[  -B+  D  B~  1 6q=  D6q* 


(10a) 
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SOCBT  Flow  at  Angle  of  Attack 

Here  the  3-D  CSCM-S  method  has  been  applied  to  predict  the  aerodynamic 
flow  field  of  a  projectile  configuration  flying  at  Mach  number  0.96,  0.98,  and  1.1  at 
4°  angle  of  attack,  and  Mach  number  0.95  at  2°  angle  of  attack.  The  predictions 
show  agreement  in  the  quantitative  and  qualitative  nature  of  the  flow  field,  including 
the  recirculatory  base  flow,  and  the  comparisons  with  experimental  data  and  other 
numerical  results.  Pressure  measurements  for  this  projectile  supported  by  a  base- 
mounted  sting  were  made  by  Kayser  and  Whiton17.  The  present  calculations  were 
done  with  a  grid  which  contained  128  longitudinal,  56  normal,  and  31  circumferential 
points.  Due  to  the  symmetry  plane  imposed  by  the  angle  of  attack  of  the  mean 
flow  and  the  axisymmetric  geometry  of  the  projectile,  only  half  of  the  domain  was 
simulated.  Figure  5  shows  the  99  by  31  grid  configuration  on  the  projectile  surface. 
The  3-D  grid  was  generated  by  rotating  the  2-D  grid  of  Figures  1  about  the  projectile 
axis  of  symmetry. 

For  a  free  stream  Mach  number  of  0.96  and  4°  angle  of  attack,  Figure  6a  shows 
the  computed  distribution  of  the  surface  pressure  coefficient  in  the  windward  and 
leeward  planes  together  with  the  square  and  circle  symbols  showing  the  experimen¬ 
tal  measurements.  The  expansions  and  recompressions  near  the  ogive-cylinder  and 
cylinder-boattail  junctions  are  well  captured  by  the  simulations.  The  agreement 
with  the  experimental  data  is  better  on  the  ogive  and  cylinder  surfaces  than  on  the 
boattail  surface. 

Figure  6b  shows  the  comparison  obtained  by  Sahu51  for  this  projectile  config¬ 
uration.  As  in  the  case  of  the  2-D  simulations,  the  upwind  method  is  in  generally 
better  qualitative  agreement  with  experiment  in  following  flow  structure  details 
than  the  central  difference  methods.  A  Mach  contour  plot  for  the  solution  of  Figure 
6a  is  shown  in  Figure  7. 

Figure  8a  shows  the  CSCM-S  computed  distribution  of  the  surface  pressure 
coefficient  in  the  leeward  plane  with  a  free  stream  Mach  number  0.95  and  2°  angle 
of  attack.  Figure  8b  shows  three  different  predictions  obtained  by  Nietubicz  et 
al52  with  different  numerical  methods  for  a  projectile  with  a  similar  configuration 
and  a  7.5°  boattail  angle.  The  CSCM  predictions  agree  better  with  the  VSYM3D 
predictions  obtained  with  the  finer  grid. 

Additional  computational  results  for  3-D  angle  of  attack  solutions  will  be  pre¬ 
sented  in  reference  36. 

The  present  initial  results  show  the  predictive  computational  aerodynamic  ca¬ 
pabilities  of  the  3-D  CSCM  method  for  projectiles  at  transonic  velocities.  As  in 
cases  run  at  higher  Mach  numbers,  the  3-D  code  predictions  are  in  quite  good 
agreement  with  results  from  the  well  validated  2-D  CSCM  codes.  This  efficient 
numerical  method  requires  about  one  hour  of  Cray-XMP  computer  time  for  each 
three-dimensional  simulation,  while  the  present  results  have  been  verified  by  using 
an  additional  hour  per  calculation.  The  method  provides  a  promising  computational 


capability  to  simulate  the  full  flow  field  around  projectiles,  including  the  separated 
base  flow. 


Various  aspects  of  computational  technique  were  explored  in  the  context  of  the 
MICOM  tactical  missile  model  base  flow  experiment15.  The  axisymmetric  problem 
was  tested  in  a  Mach  1.4  freestream.  The  model  featured  a  broad  flat  base  five 
times  the  exit  diameter  of  an  axis  centered  cold  propulsive  jet  emitted  from  the 
base.  The  jet  exit  Mach  number  was  2.7  and  the  flow  was  underexpanded  with 
normal  pressure  ratio  of  2.15  relative  to  freestream. 

We  were  particularly  interested  to  study  the  effects  of  mesh  topology  and  reso¬ 
lution  on  the  prediction  of  base  pressure,  which  was  available  from  the  experiment 
and  had  not  previously  been  successfully  simulated.  All  computational  experiments 
were  performed  with  the  same  classical  mixing  length  wake  model30. 

Figures  9a,  9b,  9c,  and  9d  show  four  meshes  that  were  studied  -  respectively 
termed  double  wraparound,  Cartesian  step  and  independent  multiple  patched,  the 
latter  without  and  with  computed  flow  structure  adaptive  grid  redistribution53.  Fig¬ 
ures  10a,  10b,  10c,  and  lOd  show  the  corresponding  Mach  contour  plots.  Comparing 
the  grid  pictures  and  the  Mach  contour  plots  one  can  clearly  see  the  correlation  be¬ 
tween  good  grid  resolution  and  flow  structure  sharpness.  The  classical  composite 
meshing  techniques,  Figures  9a  and  9b,  waste  mesh  points  in  excessive  refinement 
where,  through  the  topology,  zones  of  boundary  layer  mesh  propagate  out  into  the 
computational  domain.  The  multiple  independent  patched  mesh  strategy,  with  sin¬ 
gle  mesh  cell  overlap  to  exchange  data  by  interpolation42,  permits  better  balanced 
resolution  and  use  of  mesh  resources.  The  latter  has  small  wraparound  patches 
communicating  the  boundary  layers  around  the  sharp  corners  at  the  base. 

Finally,  Figures  11a  and  lib  show  comparisons  of  computed  base  pressure 
against  the  experimental  pressure  tap  data.  Figure  11a  shows  that  the  double 
wraparound  grid  provides  a  good  average  base  pressure  prediction,  but  a  wave 
structure  due  to  cooperative  vortex  jet  action  in  the  separated  flow  region  is  an 
artifact  of  the  calculation.  The  overall  solution  is  substantially  improved  vis  a  vis 
earlier  attempts  by  other  workers  with  central  difference  methods  and  Cartesian 
step  grids,  Figure  12  from  reference  15.  Our  results  with  Cartesian  step  grid  is  also 
somewhat  less  satisfactory  than  with  the  double  wraparound  grid. 

Figure  lib  shows  that  the  independent  multiple  grid  solution  without  adap¬ 
tation  performs  qualitatively  similar  to  the  double  wraparound  grid  in  terms  of 
the  vortex  action.  However,  the  same  topology  with  flow  structure  adaptive  grid 
point  redistribution  provides  the  additional  refinement  and,  likely  also,  alignment 
features  to  permit  a  subtle  circulation  mechanism  to  be  set  up  in  the  modified  base 
separation  vortex  pattern;  and  a  correct  flat  pressure  distribution  results.  Thus 
here  the  flow  physics  and  associated  qualitative  solution  structure  are  found  to  be 


strongly  influenced  by  grid  topology  and  resolution.  More  details  are  available  in 
references  30  and  32. 

3.  Publications 


Two  meeting  papers  in  the  area  of  algebraic  grid  generation  have  resulted,  refer¬ 
ence  33  and  the  paper  in  preparation:  Luh,  R.C.-C.  and  Lombard,  C.K.:  “FASTWO 
-  A  2-D  Interactive  Algebraic  Grid  Generator,”  AIAA-88-0516. 

Two  meeting  papers  in  the  area  of  SOCBT  flow  analyses  with  the  2-D  and  3-D 
CSCM-S  codes  have  resulted,  reference  29  and  the  paper  in  preparation:  Bardina, 
Jorge,  Lombard,  C.K.  and  Luh,  R.C.-C.:  “CSCM  Three  Dimensional  Navier-Stokes 
Computational  Aerodynamics  for  a  Projectile  Configuration  at  Transonic  Veloci¬ 
ties”,  abstract  submitted  for  AIAA  6th  Applied  Aerodynamics  Conference,  June 
1988. 

Three  meeting  papers,  references  30-32,  have  resulted  from  the  computational 
research  in  the  MICOM  propulsion  base  flow  problem. 

4.  Scientific  Personnel 

Dr.  C.K.  Lombard,  Principal  Investigator 

Dr.  Jorge  Bardina 

Dr.  Raymond  C.-C.  Luh 

Dr.  Ethiraj  Venkatapathy 
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